
#ifndef __GRID3D__H
#define __GRID3D__H    1

#include <stdio.h>


#ifdef	__cplusplus
extern "C"
{
#endif

/* define constants for different boundary conditions */
enum {
    BC_NON = 0, BC_NEU, BC_DIR, BC_SURF, BC_ZERO = 16
};
/* BC_ZERO is used to mark edges whose div-free functions are not used */

/* a grid node has three coordinates and a boundary condition */
typedef struct {
    double coords[3];
    unsigned bcond; /* Boundary Condition Code */
    unsigned on_bdry; /* bit field */
    /* =0 for interior nodes,
     * >0 for nodes on outer boundary.  If bit i is 1 node belongs to boundary 
     *    with index i.  
     */
} node_t;

/* an element edge is a line segment */
typedef struct {
    int nodes[2];
    int nf;
    int *faces;
    unsigned bcond;
#ifdef STORE_MEASURES
    double Length;
    double t[3];
#endif
} edge_t;

/* an element face is a triangle with 3 nodes, 3 edges, belongs to
    2 elements and has area and normal */
typedef struct {
    int nodes[3];
    int edges[3];
    int elems[2];
    unsigned bcond;
#ifdef STORE_MEASURES
    double Area;
    double n[3];
#endif
} face_t;

/* a finite element has four nodes, 6 edges, 4 faces and volume */
typedef struct {
    int nodes[4];
    int edges[6];
    int faces[4];
#ifdef STORE_MEASURES
    double Volume;
#endif
} elem_t;

typedef struct {
    int NN, /* Number of nodes */
            NNI, /* Number of nodes inside */
            NF, /* Number of faces */
            NFI, /* Number of faces inside */
            NE, /* Number of edges */
            NEI, /* Number of edges inside */
            NDFI, /* Number of div-free base functions inside */
            NEL; /* Number of ELements */
    char *gridname; /* The name of the grid */
    node_t *nodes; /* array of nodes */
    edge_t *edges;
    face_t *faces;
    elem_t *elems; /* array of elements */
    int *e_faces; /* array of edge faces */

    void *extended; /* different finite elements can keep data here */

    double Volume;
    double MINJAC;
    double MINQ; /* min quality of elements in grid */
    double MAXH; /* max length of edges in grid */

    struct _inverse_s *inverse; /* inverse connectivity */
} grid3d_t;


grid3d_t *GenerateGrid1(double X0, double X1, int NX,
        double Y0, double Y1, int NY,
        double Z0, double Z1, int NZ);
grid3d_t *GenerateGrid(int n);
void ReleaseGrid(grid3d_t **grid);
void LogGridStats(const grid3d_t *grid);

/* Do not call these directly unless for debugging purposes.
    They are called by GenerateGrid and LoadGrid & LoadGridBin */
void EnrichGrid(grid3d_t *grid);
void MakeFaces(grid3d_t *grid);
void SortNodes(grid3d_t *grid);
void MakeEdges(grid3d_t *grid);


typedef struct /* inverse connectivity for nodes */ {
    int NN; /* number of connected nodes */
    int NF;
    int NEL; /* number of elements containing node */
    int *nids;
    int *fids;
    int *Tids;
} node_inv_t;

/*
typedef struct
{
    int NN;
    int *nids;
} edge_inv_t;
 */
typedef struct _inverse_s {
    node_inv_t *ninv;
    /* edge_inv_t *einv; */
} inverse_t;


void MakeInverse(grid3d_t *grid);
void ReleaseInverse(grid3d_t *grid);


/****************************************************************************/
/*   The following are is g3dutils.c */
/****************************************************************************/


double GridVolume(grid3d_t *grid);

/* The following routines compute what their names refer to */
double ComputeVolume(const node_t *nodes, const int *T_nodes);
double ComputeJacobian(const node_t *nodes, const int *T_nodes);
void ComputeCentroid(const node_t *nodes, const int *T_nodes, double *pt);
double ComputeArea(const node_t *nodes, const int *f_nodes, double *n);
double ComputeLength(const node_t *nodes, const int *e_nodes, double *t);

void ComputeGlobalCoord(const node_t *nodes, const int *Tnodes,
        const double *loc, double *glo);

/* The following routines either call the corresponding routine above, or
 look up the result, depending on the STORE_MEASURES/COMPUTE_MEASURES macros */
double elem_volume(const grid3d_t *grid, const elem_t *T);
double face_area(const grid3d_t *grid, const face_t *f, double *n);
double face_vecs(const grid3d_t *grid, const face_t *f,
        double *t, double *s, double *n);
double edge_length(const grid3d_t *grid, const edge_t *e, double *t);

/* The following routine computes the quality measure of an element */
double ElemQuality(const node_t *nodes, const int *T_nodes);
/* the following routine computes the quality of elems containing node */
double NodeQuality(const grid3d_t *grid, const node_t *nodes, int nid);
/* The following routine computes the quality of the grid */
double GridQuality(const grid3d_t *grid, const node_t *nodes);
/* The following routine computes H of the grid */
double GridSize(const grid3d_t *grid, const node_t *nodes);

void UpdateMeasures(grid3d_t *grid);

/* elem_face_nodes returns a pointer to internal array of node indices,
   such that if passed to ComputeArea as f_nodes argument the returned
   unit vector n is normal to face T->faces[fn] and pointing out of T */
int *elem_face_nodes(const int *Tnodes, int fn);

/* elem_face_nodes_local returns int n[3] such that Tnodes[n[i]]==fnodes[i]
 * for i=0,1,2 */
int *elem_face_nodes_local(const int *Tnodes, const int *fnodes);


/* elem_face_area_normal returns area of face fn of element T and loads
    into NVEC the unit normal vector to this face which is outer to T.
    NOTE: fn is local index. */
double elem_face_area_normal(const grid3d_t *grid, const elem_t *T,
        int fn, double *NVEC);
/* elem_face_normal_sign returns 1 if the normal to face fn of
    element T is outer to T and -1 otherwise.
    NOTE: fn is local index. */
int elem_face_normal_sign(const grid3d_t *grid, const elem_t *T,
        int fn);

int edge_face_normal_sign(const grid3d_t *grid, const edge_t *e,
        const face_t *f);
double edge_face_coef_vect(const grid3d_t *grid, const edge_t *e,
        const face_t *f, double *V);


/* local matrices ***********************************************************/
void local_mass_p1(const grid3d_t *grid, const elem_t *T, double ***M);
void local_stif_p1(const grid3d_t *grid, const elem_t *T, double ***S);
void local_advm_p1(const grid3d_t *grid, const elem_t *T, double **AV, double ***_M);
void local_advm_D_p1(const grid3d_t *grid, const elem_t *T, double **AV, double ***_M);


void local_mass_cr(const grid3d_t *grid, const elem_t *T, double ***M);
void local_stif_cr(const grid3d_t *grid, const elem_t *T, double ***S);
void local_div_cr(const grid3d_t *grid, const elem_t *T, double ***D);
void local_div_cr_xyz(const grid3d_t *grid, const elem_t *T, double ***D);

void local_mass_cr_tsn(const grid3d_t *grid, const elem_t *T, double ***M);
void local_mass_tsdf(const grid3d_t *grid, const elem_t *T, double ***M);


void local_mass_df(const grid3d_t *grid, const elem_t *T,
        double ***M, double ***C);


/****************************************************************************/
/*   The following are is g3dio.c */
/****************************************************************************/


void print_nodes(grid3d_t *grid, const node_t *nds, char *fname);
void print_edges(grid3d_t *grid, char *fname);
void print_faces(grid3d_t *grid, char *fname);
void print_elems(grid3d_t *grid, char *fname);

void MakeGridNameFromPath(grid3d_t *grid, const char *path);

grid3d_t *LoadGridBin(const char *path);
grid3d_t *LoadGrid(const char *path);

void SaveGrid(const grid3d_t *grid, const char *filename);
void SaveGridBin(const grid3d_t *grid, const char *filename);

void save_grid_bin(const grid3d_t *grid, FILE *F);
void load_grid_bin(grid3d_t *grid, FILE *F);

#ifdef	__cplusplus
}
#endif



#endif

